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1.0 


SUMMARY 


This  report  describes  research  on  the  application  of  a  Cardinalized  Probability  Hypothesis  Density 
(CPHD)  multi-target  filter,  parameterized  in  tenns  of  specialized  geostationary  earth  orbit  (GEO) 
elements  to  estimate  the  state  of  resident  space  objects  in  the  geostationary  regime.  Justification  for 
the  use  of  GEO  elements  is  provided  through  comparisons  of  the  differential  entropy  and  Kullback- 
Leibler  divergence  for  linearized  filter  implementations  in  GEO  elements  versus  conventional  Earth 
Centered  Inertial  (ECI)  position  and  velocity  states.  Algorithms  are  demonstrated  in  a  simulation  of 
ground-based  observations  of  a  cluster  of  five  geostationary  satellites. 

The  overall  goal  of  this  work  is  to  enhance  space  situational  awareness  of  objects  in  the  high  value 
regions  of  space  populated  by  GEO-belt  communications  and  Earth-observation  satellites  and 
medium  earth  orbit  (MEO)-region  Global  Navigation  Satellite  System  (GNSS)  satellites.  By  tailoring 
estimation  methods  for  these  regions,  the  project  sought  to  improve  the  potential  for  ground-based 
telescopes  to  return  to  and  precisely  track  known  objects  and  to  detect  and  characterize  new  objects 
that  appear  in  their  field  of  view.  The  work  completed  so  far  focused  on  quantifying  the  impact  of 
the  state  representation  and  setting  up  a  CPHD  filter  to  produce  an  accurate  multi-target  solution. 


2.0  INTRODUCTION 

This  report  documents  research  performed  by  the  Colorado  Center  for  Astrodynamics  Research 
(CCAR)  at  the  University  of  Colorado  Boulder  to  enhance  space  situational  awareness  (SSA)  through 
the  use  of  multi-target  filtering  methods.  The  work  focused  on  the  application  of  Finite  Set  Statistics 
(FISST)  and  a  specialized  state  representation  to  estimate  the  orbits  of  multiple  space  objects  in  the 
GEO  orbital  regime.  We  had  planned  to  also  look  at  MEO;  however,  in  the  available  time  we  were 
only  able  to  complete  the  investigation  of  the  GEO  case. 

2.1  GEO  Elements 

It  is  well  known  that  the  performance  of  linearized  estimation  methods  such  as  the  Extended  Kalman 
Filter  (EKF)  can  be  degraded  when  the  true  system  dynamics  or  observation  dependence  on  the 
system  states  are  not  well-approximated  by  linear  models  and  many  researchers  have  explored 
parameterizations  that  extend  the  applicability  of  linear  models  [1].  If  these  issues  are  not  properly 
addressed,  use  of  linearized  filters  can  produce  inaccurate  state  estimates  and  error  covariance 
predictions  when  the  nonlinearities  produce  probability  distributions  of  the  observation  or  state  errors 
that  are  not  Gaussian.  Weisman  et  al.  [2]  looked  specifically  at  a  rigorous  technique  for  mapping 
uncertainties  in  nonlinear  systems.  Many  references  also  provide  mapping  functions  to  facilitate 
conversion  of  states  and  covariance  matrixes  from  one  representation  to  another;  Vallado  [3]  provides 
a  convenient  summary  of  these. 
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Significant  research  has  also  been  performed  to  develop  estimators  that  are  more  effective  in  dealing 
with  nonlinearities  and  non-Gaussian  measurement  errors.  The  most  prevalent  is  probably  the 
Unscented  Kalman  Filter  (UKF)  [4,  5].  Despite  the  existence  of  methods  like  the  UKF,  the  possibility 
of  performing  estimation  using  a  simple  linearized  filter  is  still  of  interest  if  it  can  be  shown  to  produce 
reliable,  realistic  results  with  lower  computational  complexity.  Thus,  the  first  part  of  this  work 
investigated  the  benefit  of  using  a  specialized  state  representation,  tenned  GEO  elements  [6,  7],  which 
reduce  both  the  nonlinearity  of  the  state  propagation  and  the  relationship  between  the  state  and  angles- 
based  observations  of  the  object  from  a  ground  station.  These  are  compared  with  conventional  ECI 
position  and  velocity  state  vectors  [8]. 

2.2  Measures  of  Performance 

Two  measures  were  used  to  compare  the  performance  of  a  GEO  element  state  vector  representation 
to  a  conventional  ECI  state  vector  representation  in  EKF  and  UKF  filter  implementations.  The 
Kullback-Leibler  (KL)  divergence  [9]  provides  a  means  to  compare  two  Probability  Density 
Functions  (PDF),  and  identify  when  information  is  gained  in  an  estimator.  The  Differential  Entropy 
(DE)  metric  has  also  been  used  in  particular  to  determine  when  an  error  distribution  becomes  non- 
Gaussian  [10]. 

Both  measures  show  that  the  GEO  representation  can  be  accurately  propagated  and  used  to 
incorporate  measurements  using  a  simple  EKF  implementation  without  a  meaningful  loss  of 
information. 

2.3  FISST  Methods 

Using  the  GEO  element  representation,  a  multi-target  CPHD  filter  [11,  12,  13]  was  implemented, 
with  the  goal  of  tracking  multiple  objects  observed  by  a  ground-based  sensor  in  the  presence  of 
significant  clutter.  The  perfonnance  of  the  CPHD  was  tested  using  a  simulation  of  a  particular 
commercial  GEO  satellite  cluster  of  five  satellites  sharing  a  single  orbital  slot.  A  future  effort  will 
apply  the  techniques  to  experimental  observations  of  a  larger  number  of  targets. 

3.0  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

This  section  provides  the  background  models  and  methods  used  for  the  study  related  to  the  GEO 
representation,  the  performance  measures,  and  a  brief  introduction  to  the  CPHD  multi-target  tracking 
approach. 


3.1  State  Parameterization-GEO  Elements 

A  key  objective  of  this  work  was  to  utilize  state  parameterizations  that  optimize  linear  propagation. 
Toward  this  end,  the  following  set  of  non-dimensional  orbital  elements,  optimized  for  GEO,  were 
chosen  [6,  7]: 


a  —  +  n  +  v)  —  fi(t)  (i) 
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a  —  A 

A  a  -  - 

A 

(2) 

ex  —  e  cos(n)  +  CL) 
ey  —  e  sin(m  +  fl) 

(3) 

(4) 

Qi  =  tan  Q  sin  (13) 

(5) 

Q2  =  tan  Q  cos (13) 

(6) 

Where  the  GEO  elements  are  defined  in  terms  of  the  Keplerian  elements:  semimajor  axis  a, 
eccentricity  e,  inclination  i,  right  ascension  of  the  ascending  node  CL,  argument  of  perigee  <n,  and  true 
anomaly  v.  The  GEO  elements  are:  Earth- fixed  sub-satellite  longitude  A,  non-dimensional 
longitudinal  drift  rate  A  a,  eccentricity  vector  components  (ex,  eyJ,  and  the  equinoctial  elements 
(Qi’Qz)-  The  longitudinal  drift  rate  is  nonnalized  by  a  nominal  GEO  semimajor  axis:  A  — 
42,164.2  km.  Finally,  the  longitude  includes  the  Greenwich  sidereal  angle  6(t),  which  is  a  function 
of  the  average  Earth  rotation  rate  u)E: 

G(t)  =  G0  +  (oE(t  —  t0)  (7) 

Computation  of  the  osculating  GEO  elements  from  known  ECI  position  and  velocity  vectors  can  be 
completed  by  first  forming  the  classical  Keplerian  elements  and  then  using  (l)-(7)  above  to  find  the 
new  element  set.  Vallado  ([8]  RV2COE,  p.  1 13)  lays  out  the  following  procedure  to  get  the  Keplerian 
elements: 


h  —  r  x  r 

(8) 

h  =  \h\ 

(9) 

n  —  K  x  h 

(10) 

2  “  r)  f  “  (f '  ^ 

(11) 

k 

e  =  \e\ 

(12) 

f2  lA 
£  =  --- 
2  r 

(13) 

/* 

a  —  —  — 

2£ 

(14) 

h7 

cos  (i)  =  — 
n 

(15) 
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co s(I2)  — if  ny  <  0,  then  fl  =  360°  —  fi  (16) 

cos(u>)  =  if  ez  <  0,  then  co  =  360°  —  co  (17) 

cos(v)  =  if  f  ■  f  <  0,  then  v  —  360°  —  v  (18) 

Note  that  other  equivalent  definitions  are  also  available  in  the  literature. 

3.2  State  Estimate  Covariance 

Under  the  assumption  of  Gaussian  error  distributions,  the  state  error  covariance  matrix  provides  a 
means  to  compare  the  perfonnance  of  different  approaches.  This  report  explores  the  characteristics 
of  the  covariance  matrix  for  two  different  state  representations.  To  compare  them  an  accurate 
transformation  is  required. 

In  its  most  simple  form,  the  covariance  matrix  is  diagonal  with  the  squares  of  the  standard  deviation 
of  the  state  component  errors  along  the  diagonal.  For  the  ECI  state  representation,  it  might  be  given 

by, 


P ECI  ~ 


\  0 

0  Gy 

0  0 

0  0 

0  0 

.  o  0 


0  0 

0  0 

erf  0 

0  4 

0  0 

0  0 


0  0  ■ 
0  0 
0  0 
0  0 
of  0 

0  of. 


(19) 


In  general  there  are  also  off-diagonal  terms  indicating  correlations  between  the  states.  A  similar 
covariance  matrix  can  be  formed  for  the  GEO  elements.  In  our  work,  two  methods  were  utilized  to 
transform  between  state  representations:  the  unscented  transform  and  the  matrizant.  The  unscented 
transform  requires  generating  a  set  of  sigma  points  to  represent  the  original  state  distribution, 
transforming  each  of  these  points  via  a  full  nonlinear  mapping  into  the  desired  coordinates,  then 
recalculating  the  covariance  matrix  in  the  new  coordinates  by  differencing  all  the  transformed  sigma 
points  with  the  nominal  state.  The  matrizant  or  similarity  transform  is  typically  less  computationally 
intensive,  but  does  require  the  calculation  of  partials  of  one  state  representation  with  respect  to  the 
other.  As  mentioned  above,  analytical  partials  are  available  for  a  number  of  common  transformations 
in  [3],  The  two  methods  are: 


PGE0  —  Unscented  J'ransform(PECl ) 


(20) 


Pgeo 


dXGE0  p  dXGE0 
dXEGI  dXEGI 


(21) 


4 

Approved  for  public  release;  distribution  is  unlimited. 


3.3 


GEO  Element  Propagation  Equations 


A  key  advantage  expected  for  use  of  the  GEO  elements  is  the  linearity  and  slower  time  rate  of  change 
of  the  elements.  Nonlinear  GEO  element  variational  equations  in  terms  of  the  body-fixed  radial  (ar), 
tangential  (ae),  and  out-of-plane  (ah)  disturbing  forces  were  derived  and  presented  in  Refs  [6,  7],  as 
follows: 


h  r 

2  =  +  -  [Q2  sin(s)  -  Q1  cos( s)]ah  -  uE 

rz  h 


2 (5u  +  l)2  r,  .  p 

Sa  — - — - \(ex  sin(s)  —  ey  cos(s)Jar  +  — 


hA 


r  (V  f  V\ 

—  j—  sin(s)  ar  +  ex  +  (^1  +  — J  cos(s)  ae  +  ey  [Qi  cos(s)  - 


r—p  f  V\ 

| — cos(s)  ar  +  ey  +  (^1  +  -J  sin(s)  ae  —  ex[Q1  cos(s)  — 


Qi  =  2^(!  +  Ql  +  Ql)  sin (s)  ah 


=  2h  C1  +  ^  +  ^  C°S^  Uh 


(22) 

a-e] 

(23) 

■  Q2  sin(s)]ah} 

(24) 

■  Q2  sin(s)]ah} 

(25) 

(26) 

(27) 

Comparable  variational  equations  for  ECI  position  and  velocity  vector  components  are  readily  found 
in  Vallado  ( [3],  Numerical  Integration  p.  591). 


3.4  Filter  Update  Equations 

Propagation  of  the  state  estimates  is  perfonned  by  numerical  integration  of  the  equations  of  motion 
given  above  for  the  GEO  states  and  in  [8]  for  the  ECI  states.  For  this  study,  the  integration  was 
performed  using  built-in  integration  methods  in  MATLAB.  The  propagation  of  the  covariance  is 
performed  in  two  different  ways  -  assuming  linear  updates,  as  in  the  EKF,  or  via  nonlinear  methods 
via  sigma  points  UKF.  The  linear  EKF  time  update  for  the  covariance  is: 

P(t)  =  «D(tkf  tk_i)P(t  -  l)d>(tk,tk_1)7’  +  Q  (28) 

where  <t)(tk,  tk-i)  is  the  state  transition  matrix. 

The  measurement  update  for  the  EKF  uses  the  measurement  Jacobian  H  with  measurement  noise 
variance  R : 


S  =  HP~Ht  +  R  (29) 
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To  calculate  the  Kalman  gain  K: 


K  =  p-tfS-1  (30) 

Giving  the  linear  update  equation  for  the  covariance: 

Pl+inear  =  06x6-KH)P-  (31) 

The  superscript  (-)  denotes  the  pre-update  covariance  and  (+)  is  the  post-update  covariance. 

Following  [4,  5],  nonlinear  propagation  in  the  UKF  time  updates  requires  integrating  the  equations 
of  motion  for  13  sigma  points,  giving  the  mean  state  and  covariance  as: 

N 

Z(t)  =YjWMiXi(t)  (32) 

i= 1 


N 

p(t) = ^  waiXiit)  -  mmct)  -  mY  w 

i= 1 

Where  WM  and  Wc  are  weights  for  the  mean  and  covariance  calculations  respectively.  The  sigma 
points  Xt  (t)  are  propagated  using  a  numerical  integrator.  And  X(t)  is  the  mean  of  the  sigma  points. 
Measurement  updates  also  require  the  generation  of  sigma  points.  Computed  measurements  are 
generated  through  a  measurement  function  h(-)  with  error  £{. 

Z,  =  h{Xd  +  et  (34) 

The  generation  of  sigma  points  includes  the  generation  of  weights  for  each  sigma  point  Wxi,  allowing 
us  to  add  up  the  computed  measurements  and  get  a  mean  value  Z : 

N 

Z  =  WmZi  (35) 
i= 1 

Calculating  two  intermediate  variables: 

N 

5  =  ^VFa(Z,-Z)(Zi-Z)r  (36) 

i= 1 
N 

G  =  ^  Wci  (Xi  -  X)  C Zf  -  zy  (37) 

i= 1 

Allows  the  update  equation  under  the  UKF  paradigm: 
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Pnonlin  =  P~  ~  GS~1GT  (38) 

3.5  Process  Noise  Discussion 

Much  effort  was  invested  in  being  certain  that  process  noise  was  added  equivalently  in  the  two 
coordinate  frames.  For  the  simulations  in  this  work  the  initial  covariance  and  process  noise  matrix 
Q  are  specified  in  Radial,  In-track,  Cross-track  (RIC)  coordinates.  They  are  then  rotated  into  ECI. 

P  —  rRlc  p  rRicT  (39t 

rECI  —  lEC1  rRlC  lEC1 

For  the  filters  running  in  GEO  elements,  the  matrices  undergo  a  similarity  transfonnation  to  GEO 
element  space. 


Pgeo  — 


dX, 


GEO 


dX, 


GEO 


dX 


ECI 


dX, 


ECI 


Alternatively  the  conversion  may  be  done  with  an  unscented  transfonn. 


(40) 


3.6  Differential  Entropy 

In  Reference  [10],  DeMars,  et  ah,  introduced  DE  as  a  means  to  adapt  a  Gaussian  Mixture  Model 
(GMM)  to  correctly  represent  non-Gaussian  state  uncertainties  that  arise  from  nonlinearity  in  state 
propagation. 

The  general  equation  for  DE  is  defined  as: 

H(x)  =  -  J p(x)  log (p(x))  dx  -  E[-  log(p(x))]  (41) 

where  p(x)  is  the  probability  density  of  the  state,  and  the  expected  value  is  taken  with  respect  to 
p(x).  Specifically  applied  to  a  Gaussian  probability  density  function  (PDF)  this  becomes  [10]: 

1 

H(x)  = -\og\2neP\  (42) 
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Where  the  mathematical  constant  is  e  =  2.71828  and  P  is  the  covariance  matrix  describing  the 
Gaussian  distribution.  The  entropy  metric  describes  the  average  amount  of  information  content 
associated  with  a  random  variable. 


By  comparing  the  differential  entropy  of  a  filter  before  and  after  a  time  propagation,  or  before  and 
after  a  measurement  update,  one  can  assess  the  loss  or  gain  of  information.  Furthennore,  by 
comparing  the  differential  entropy  for  a  linearized  EKF  filter  to  the  UKF  filter,  the  impact  of  the 
linearization  can  be  assessed.  We  use  this  approach  to  compare  filters  implemented  in  GEO 
elements  to  a  Cartesian  state  representation.  One  concern  is  that  if  the  distribution  is  not,  in  fact, 
Gaussian  or  the  filter  computed  covariance  matrix  is  not  accurate,  using  Equation  40  in  this  manner 
might  not  be  justified. 


3.7  Kullback-Leibler  Divergence 

We  also  investigated  the  use  of  KL  divergence  as  a  measure  of  degradation  in  the  uncertainty  due 
to  unmodeled  non-linearities  [9].  Also  known  as  the  information  divergence,  KL  is  a  measure  of 
the  difference  between  two  probability  distributions. 

To  assess  the  linearity  of  a  proposed  state  representation,  the  KL  divergence  can  be  used  to  compare 
the  results  using  linear  approximations  in  an  EKF  (PEKF),  with  non-linear  implementation  in  a  UKF 
(PlJKFy  if  qlc  siaie  ts  we]|  represented  by  the  linear  methods,  one  would  expect  to  see  little  difference 
between  the  results.  The  KL  divergence  can  also  compare  the  information  change  due  to  a 
measurement  update. 


The  KL  divergence  for  a  Gaussian  PDF  may  be  computed  as  [9]: 


Dkl  ( Pukf\\Pekf ) 


(43) 


The  KL  divergence  in  this  form  comprises  the  sum  of  three  terms,  highlighted  in  the  bracket  in  the 
equation  above.  The  first  term  is  referred  to  here  as  the  trace  effect.  The  additional  variable  K  is 
simply  the  dimension  of  the  state  vector,  always  six  for  this  work.  In  the  case  where  the  two 
covariance  matrixes  are  identical,  the  trace  of  the  product  of  the  covariance  matrix  and  its  inverse 
will  be  equal  to  its  dimension,  so  by  subtracting  this  value,  trace  effect  term  goes  to  zero.  The  second 
tenn  compares  the  actual  value  of  the  state  errors  against  the  expected  covariance.  Finally,  the  last 
tenn  can  be  viewed  as  a  volume  metric,  comparing  the  volumes  of  the  two  error  hyper-ellipsoids.  If 
the  covariance  matrixes  are  identical,  this  term  will  also  go  to  zero. 
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3.8 


CPHD  Description 


The  end  goal  of  our  work  is  to  use  GEO  elements  to  better  track  multiple  targets  for  GEO  space 
situational  awareness.  The  application  of  FISST  based  techniques  to  SSA  has  been  receiving 
considerable  attention  in  the  literature  because  of  its  advantages  in  terms  of  removing  the  requirement 
for  association  of  observations  with  specific  objects  and  establishing  a  rigorous  framework  for  large- 
scale  observation  problems  with  multiple  targets  and  sensors. 

Of  particular  interest  is  the  cardinalized  probability  hypothesis  density  CPHD  filter,  a  powerful 
approach  for  multi-target  tracking  offering  that  in  addition  to  the  overall  FISST  benefits,  also  offers 
a  stable  estimate  of  the  number  of  objects  that  is  useful  in  its  own  right,  and  provides  useful  bounds 
for  large-scale  observation  problems  like  SSA.  The  challenge  with  the  CPHD  is  that  in  general  it  is 
not  tractable;  however,  reference  12  derived  an  analytical  implementation  of  the  CPHD  filter 
applicable  to  linear  (and  nearly  linear)  Gaussian  models.  Thus,  if  the  GEO  state  representation 
proposed  above,  does  in  fact  meet  the  requirements  for  linearity,  the  Gaussian  Mixture  (GM)  CPHD 
promises  to  be  an  excellent  approach  for  SSA  of  geostationary  objects. 

Details  of  the  GM  CPHD  derivation  and  implementation  are  provided  in  [1 1,  12].  We  followed  the 
procedures  from  that  reference,  with  additional  insights  and  advice  from  B.  Jones  and  S.  Gehly.  The 
foundation  of  the  filter  is  the  use  of  FISST  and  Gaussian  mixtures  to  efficiently  represent  the  state 
and  to  consider  every  possible  observation  association  hypothesis.  At  each  measurement  update,  the 
weight  associated  with  a  particular  hypothesis  increases  if  an  observation  can  be  matched  with  a 
predicted  target  state,  otherwise  the  weight  decreases.  Once  a  weight  decreases  beyond  some 
threshold,  the  hypothesis  is  discarded.  New  Gaussian  components  are  added  to  the  state  via  different 
models:  birth  of  targets  at  a  location,  spawning  of  targets  from  a  surviving  target,  and  splitting  of  a 
GM  if  the  distribution  is  found  to  no  longer  be  Gaussian.  We  have  not  yet  included  birth  and 
spawning  models  in  our  work. 

The  CPHD  recursion  comprises  a  prediction  (time  propagation)  and  a  measurement  update  of  the 
intensity  function  'tr(x)  and  the  cardinality  distribution  p(n).  In  the  GM  CPHD  both  are 
constructed  from  appropriately  weighted  single  target  models  and  estimates.  Details  are  given  in 
[11,12], 

It  should  be  noted  that  while  the  CPHD  can  estimate  target  states  based  on  dynamical  models  and 
measurement  updates,  it  does  not  explicitly  ‘track’  the  targets  in  the  sense  that  there  are  no  inherent 
labels  to  monitor  the  evolution  in  time  of  a  specific  target.  Newer  implementations  that  address  this 
issue  are  the  labeled  CPHD  filter  [14]  or  labeled  Bernoulli  filter  [15]  designed  specifically  to  track 
targets. 

To  evaluate  the  state  estimation  accuracy  in  our  preliminary  work,  we  manually  assigned  each 
estimated  object  to  the  closest  true  object,  without  duplication.  In  [13]  the  optimal  sub-pattern 
assignment  (OSPA)  metric  [16]  is  used  to  compare  the  estimated  state  to  truth.  OSPA  assigns  the 
largest  possible  subset  of  true  and  estimated  states,  then  adds  an  offset  for  all  points  that  are 
unassigned.  This  approach  will  be  implemented  in  future  extensions  of  this  work. 
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4.0 


RESULTS  AND  DISCUSSION 


4.1  Linear  Propagation  Comparison 

A  representative  initial  covariance  for  a  geostationary  satellite  was  propagated  with  nonlinear  UKF 
and  linear  EKF  time  updates  in  both  ECI  and  GEO  element  space.  The  initial  covariance  was 
specified  as  diagonal  in  ECI  and  an  unscented  transform  was  used  to  obtain  the  initial  covariance  in 
GEO  elements.  The  initial  covariance  was  set  in  ECI  with  aP  —  1  km  in  each  position  coordinate  and 
ov  =  1x10 ~4  km/s  in  each  velocity  component. 

Included  in  the  state  propagation  were  perturbing  accelerations  due  to  an  8  x  8  Earth  gravity  field, 
Solar  and  Lunar  gravity,  and  Solar  Radiation  Pressure  (SRP).  For  the  UKF,  the  sigma  points  were 
propagated  using  a  Runge-Kutta  4th  Order  Integration  Scheme  (RK4)  integrator  with  fixed  time  steps 
of  10  minutes.  The  sigma  points  were  not  re-initialized  after  each  time  step,  that  is,  one  set  of  sigma 
points  was  propagated  for  the  entire  48  hour  simulation. 

Random  disturbing  accelerations  were  modeled  using  process  noise,  defined  in  the  orbit  local  or  RIC 
frame  integrated,  and  rotated  or  transformed  (using  the  unscented  transform)  in  to  the  ECI  and  GEO 
representations. 

The  differential  entropy  was  calculated  for  the  ECI  covariance  from  both  propagations  methods, 
and  the  difference  in  DE  between  the  sigma  point  UKF  and  state  transition  matrix  EKF 
propagations,  {H{P^ci  (0)  —  W(P|cf  (0)  ),  is  plotted  in  Figure  1.  Similarly,  Figure  2  shows  the 
same  differential  entropy  metric  for  propagation  in  the  GEO  element  space,  {H  {P'geo  (0)  — 


Figure  1.  Difference  Between  Differential  Entropy  for  ECI  Propagation  Using  UKF  and  EKF 

Methods 


10 

Approved  for  public  release;  distribution  is  unlimited. 


Time  Step,  hrs 


Figure  2.  Difference  Between  Differential  Entropy  for  GEO  Element  Propagation  Using  UKF 

and  EKF  Methods 

Figure  1  shows  that  with  an  ECI  state  representation  the  DE  difference  between  a  linearized  EKF  and 
nonlinear  UKF  based  propagation  is  very  small  for  about  10  hours,  then  increases  steeply  especially 
in  the  time  between  30  and  40  hours  to  a  non-dimensional  value  just  under  2.  The  fact  that  the 
difference  is  positive  throughout,  indicates  that  the  linear  version  is  underestimating  the  uncertainty 
as  compared  to  the  nonlinear  version.  On  the  other  hand,  the  DE  difference  for  the  GEO  state 
representation  in  Figure  2  is  negligibly  small  (less  than  0.0004)  throughout  the  48  hours.  This 
outcome  is  consistent  with  prior  research  that  showed  that  equinoctial  elements,  which  are  very 
similar  to  GEO  elements,  remain  close  to  linear  for  substantially  longer  propagation  times  than 
Cartesian  coordinates. 

As  expected,  these  results  show  that  we  can  confidently  propagate  the  covariance  of  GEO  satellites 
states  with  a  state  transition  matrix  in  the  GEO  element  space,  without  losing  accuracy.  This 
simplification  reduces  the  computational  burden  of  numerically  propagating  sigma  points  in  a  UKF 
implementation  and  offers  the  potential  for  applicability  of  the  Vo  (GM-CPFID)  implementation.  The 
same  is  not  necessarily  true  for  ECI  states.  However,  it  can  be  argued  that  as  long  as  the  propagation 
interval  is  less  than  10  hours  (for  geostationary  orbits),  the  results  above  show  that  the  state  transition 
matrix  propagation  for  ECI  is  also  valid. 

4.2  Single  Measurement  Update  Comparison 

To  investigate  the  impact  of  nonlinearities  in  the  measurement  updates,  we  can  look  at  the  differential 
entropy  metric  before  and  after  a  measurement  update  is  performed  via  both  EKF  and  UKF 
formulations.  The  measurements  are  assumed  to  be  azimuth  and  elevation  angles  from  a  ground 
station.  For  this  simple  test,  and  the  analyses  in  later  sections,  the  ground  station  is  located  in 
Hartebeesthoek,  South  Africa.  The  same  initial  covariance  in  ECI  Equation  19  and  in  GEO  elements 
via  the  unscented  transform  are  used.  The  nonlinear  update  is  performed  via  sigma  point  generation 
(Aj,  fj)  and  subsequent  calculation  of  the  measurement  update.  Time  is  ignored  below  because  only 
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one  time  step  is  considered  here.  Results  of  this  single  point  test  are  tabulated  in  Table  1  below.  The 
initial  differential  entropies  in  the  first  column  are  H0  from  Equation  40.  The  second  column  shows 
the  differential  entropy  after  a  linear  EKF  measurement  update.  In  the  third  column  are  the  results 
from  the  nonlinear  UKF  measurement  update.  The  fourth  column  gives  the  change  in  DE  across  the 
measurement,  computed  as  the  difference  between  the  DE  from  the  post-measurement  UKF  (column 
3)  and  the  DE  from  the  a  priori  covariance  (column  1).  The  first  row  of  the  table  gives  the  ECI  results 
and  the  second  row  gives  the  GEO  element  results.  Notice  that  the  initial  DE  values  are  different 
between  these  rows  because  the  units  of  the  covariance  matrixes  in  the  two  representations  are 
different.  The  third  row  shows  the  DE  based  on  the  covariance  matrix  from  the  GEO  filter  (second 
row),  converted  to  ECI  using  an  unscented  transformation.  The  fact  that  the  first  column  results  in 
rows  1  and  3  are  the  same,  indicates  that  this  transformation  has  been  done  consistently. 

Comparing  columns  2  and  3  we  see  that  for  this  particular  initial  covariance  matrix  and  measurement 
set,  there  is  no  meaningful  difference  between  the  decrease  in  differential  entropy  (or  information 
gained)  produced  by  the  EKF  and  UKF  measurement  updates  for  either  the  ECI  or  GEO  filter.  This 
indicates  that  the  measurement  is  quite  linear  with  respect  to  both  state  representations. 

Column  4  of  Table  1  shows  that  the  same  information  gain  is  achieved  by  the  measurement  in  both 
ECI  and  GEO  when  a  measurement  update  is  performed.  Comparison  of  the  column  4  values  in  rows 
2  and  3  shows  that  the  perfonnance  metrics  are  preserved  when  the  results  are  transformed  from  GEO 
to  ECI. 

While  interesting,  it  is  important  to  note  that  Table  1  represents  a  single  test  case.  Additional  results 
and  analysis  are  required  to  fully  quantify  in  the  more  general  case,  how  much  information  is 
gained/lost  in  making  measurement  updates  in  GEO  versus  ECI. 

Table  1.  Differential  Entropy  Before  and  After  a  Measurement  Update 


H(P~ ) 

H(Pnonliv) 

««„„„■„)  -  H(P-) 

ECI 

-19.117 

-23.285 

-23.285 

-4.167 

GEO 

-55.129 

-59.296 

-59.296 

-4.167 

GEO  >  ECI 

-19.117 

-23.285 

-23.285 

-4.167 
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4.3 


Full  Simulation  Entropy  Comparison 


We  next  set  out  to  show  the  effects  of  different  estimators  and  state  representations  in  the  4-day 
simulation  by  saving  the  estimates  and  covariance  matrix  at  every  time  step  before  and  after  a 
measurement  update.  For  consistency  all  four  filters  were  run  on  the  same  simulated  observation 
data  with  process  noise  values  set  in  RIC  coordinates  and  transformed  to  ECI  and  GEO  as  required. 

Figure  3  presents  the  change  in  DE  over  each  of  the  measurements  for  the  ECI  (Red)  and  GEO  (Blue) 
implementations  run  with  EKF  and  UKF  including  process  noise.  For  each  case  the  top  sections 
show  the  EKF  (+)  and  UKF  (o)  DE  values,  and  the  bottom  sections  show  the  difference  between 
them.  To  highlight  the  measurement  epochs,  the  x-axis  in  these  plots  is  the  observation  number  in 
the  simulation,  not  the  elapsed  time.  In  the  top  plot  the  four  results  are  indistinguishable.  In  the 
bottom  plot  one  can  observe  very  small  discrepancies  in  DE  between  EKF  and  UKF  for  both  the  ECI 
and  GEO  parameterizations,  all  less  than  5E-4. 

Figure  4  shows  the  difference  in  DE  between  EKF  and  UKF  for  each  time  step  in  the  simulation 
compared  to  the  starting  value.  Subtracting  the  first  value  removes  the  unit  dependence  of  the 
covariance  matrix  on  the  DE  result.  Again  we  show  ECI  (Red)  and  GEO  (Blue)  runs  with  EKF  (+) 
and  UKF  (o).  In  the  top  plot  all  four  filter  results  are  indistinguishable.  What  we  see,  as  expected,  is 
a  decrease  in  the  DE  as  the  filters  incorporate  measurements  and  improve  their  estimates.  In  the 
bottom  plot,  shown  in  Blue,  is  the  extremely  small  difference  between  the  DE’s  computed  from  the 
GEO  EKF)  and  GEO  UKF  covariance  matrixes,  indicating  that  the  linearized  GEO  implementation 
is  valid.  The  results  in  Red  for  the  ECI  estimator  are  also  quite  small,  showing  only  a  slightly  negative 
offset  of  about  -0.015.  The  negative  sign  indicates  that  the  EKF  is  very  slightly  underestimating  the 
covariance  compared  to  the  UKF.  During  the  time  propagations  there  is  virtually  no  change  in  the 
DE  values  for  either  filter. 

The  result  of  this  study  show  that  both  ECI  and  GEO  are  suitable  for  measurement  updates  and  time 
propagation  of  states  and  covariances. 
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Observation  Number 


Figure  3.  Change  in  Differential  Entropy  Across  Each  Measurement  Update 
for  ECI  (Red)  and  GEO  (Blue)  Implementations  of  EKF  (+)  and  UKF  (o) 
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Figure  4.  Change  in  Differential  Entropy  Overtime  for  ECI  (Red)  and 
GEO  (Blue)  implementations  of  EKF  and  UKF 
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4.4 


Kullback-Leibler  Comparison 


Next  we  looked  at  the  results  using  the  KL  divergence  criteria.  Figures  5  and  6  show  the  comparison 
of  the  KL  measure  based  on  the  covariance  and  state  errors  from  the  UKF  to  the  covariance  and  state 
errors  from  the  EKF  for  the  ECI  and  GEO  state  filters.  Both  figures  show  very  little  nonlinearity, 
with  the  ECI  results  all  below  0.06  and  the  GEO  results  all  less  than  0.0012.  The  small  differences 
in  the  state  estimates  between  EKF  and  UKF  are  shown  in  Figures  7  and  8  for  the  ECI  and  GEO 
cases,  respectively. 


Figure  5.  KL  Measure  Comparing  UKF  Covariance  to  EKF  Covariance 

for  ECI  State  Filters 
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Figure  6.  KL  Measure  Comparing  UKF  Covariance  to  EKF  Covariance 

for  GEO  State  Filters 
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Figure  7.  Difference  in  the  ECI  State  Estimates  Between  EKF  and  UKF  (km  and  km/s) 
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Figure  8.  Difference  in  the  GEO  State  Estimates  Between  EKF  and  UKF  (Unit  Less) 


4.5  CPHD  Verification 

To  verify  that  the  GEO  element  EKF  version  of  the  CPHD  was  implemented  correctly,  a  single 
target  was  filtered  in  the  EKF  CPHD  and  compared  to  a  single  target  filtered  with  an  established 
filtering  software  -  Orbit  Determination  Toolbox  (ODTBX1).  For  this  comparison,  the  same 
dynamics  and  measurement  models  were  used  in  each  filter,  but  the  filter  implementations  were 
set  up  independently.  The  results  shown  in  Figures  9  and  10,  illustrate  that  the  two  are  consistent; 
thus  we  are  confident  that  the  CPHD  filter  code  has  been  implemented  correctly. 


1  Orbit  Determination  Toolbox  (ODTBX),  open  source  MATLAB  software  created  at  NASA  GSFC,  which  is  available 
online  at  http://sourceforge.net/projects/odtbx/ 
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Figure  9.  Comparison  of  ODTBX  EKF  Filtering  Results  (Left)  Versus  EKF  CPHD  (Right)  in 

GEO  Elements,  First  3  Elements 


Figure  10.  Comparison  of  ODTBX  EKF  Filtering  Results  (Left)  Versus  EKF  CPHD  (Right)  in 

GEO  Elements,  Last  3  Elements 
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4.6 


Tracking  Multiple  Simulated  GEO  Satellites 


To  put  the  multi-target  CPHD  filter  in  GEO  elements  to  the  test,  we  considered  the  observation  and 
tracking  of  an  example  set  of  GEO  satellites  sharing  the  same  orbital  slot.  The  Astra  communication 
satellites  2A,  2C,  2E,  2F,  and  2G  were  all  collocated  in  the  same  GEO  slot  as  of  June  25,  2015.  The 
Two  Line  Elements  (TLEs)  for  these  five  satellites  were  downloaded  from  CelesTrak2  and  used  to 
seed  the  initial  conditions  for  the  simulations  presented.  Initial  conditions  in  ECI  were  generated  from 
the  TLEs,  then  propagated  using  Java  Astrodynamics  Toolkit  (JAT)3  to  create  a  truth  ephemeris.  The 
forces  included  a  2 1  x  2 1  Earth  gravity  field,  Solar  gravity,  Lunar  gravity,  and  solar  radiation  pressure. 
Two  days  of  truth  data  were  generated  and  saved  in  both  ECI,  earth  centered,  earth  fixed  (ECEF), 
and  GEO  element  coordinates  (ECEF  to  facilitate  generation  of  measurements).  The  integration  was 
performed  with  RK4  formulation  and  a  five  minute  fixed  step  integration. 

Simulated  measurements  of  azimuth  and  elevation  were  generated  at  5  minute  intervals  over  the  two- 
day  simulation  from  a  Ground  Station  in  Hartebeesthoek,  South  Africa.  Measurement  errors  of  1  arc 
second  (1 -sigma)  Gaussian  white  noise  were  added  to  the  true  measurements.  In  addition  to  the  5  sets 
of  measurements  corresponding  to  the  actual  satellites,  every  measurement  update  also  included  false 
measurements  or  clutter.  Clutter  was  generated  with  a  uniform  distribution  in  azimuth  and  elevation, 
spanning  the  sensor  field  of  view  close  to  the  Astra  satellites.  The  number  of  clutter  measurements 
at  each  observation  epoch  was  determined  by  a  Poisson  distribution.  A  full  day  of  measurements  for 
all  five  satellites,  shown  in  Blue,  and  the  clutter,  shown  in  Red,  can  be  seen  in  Figure  11.  Note  that  if 
the  actual  measurements  were  not  plotted  in  a  different  color,  it  would  be  nearly  impossible  to 
distinguish  them  in  the  midst  of  the  clutter. 


Azimuth 


Time,  hrs 

Figure  11.  Measurements  of  Astra  Satellites  (Blue)  and  Clutter  (Red)  Input  to  the  Filter 


2  Website  for  downloads:  www.celestrack.com 

3  Java  Astrodynamics  Toolkit  (JAT),  open  source  Java  software  available  online  at  http://jat.sourceforge.net/ 
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Within  the  CPHD,  the  GEO  states  are  propagated  by  numerical  integration  of  the  full  nonlinear 
differential  equations;  the  covariance  is  updated  via  the  state  transition  matrix,  and  the  measurement 
updates  are  based  on  the  EKF  implementation.  To  evaluate  the  accuracy  of  the  estimates,  a  simple 
target  association  algorithm  is  used  to  match  an  estimated  target  with  a  true  target.  The  five  sets  of 
plots  in  Figure  12  show  the  errors  in  the  GEO  elements  between  the  closest  estimated  target  to  each 
simulated  Astra  satellite. 

Overall  our  initial  implementation  of  the  GEO  CPHD  filter  tracks  the  targets  reasonably  well,  given 
the  cluttered,  noisy  measurements  provided.  There  is  also  clearly  room  for  improvement.  The  a 
priori  covariance  was  set  somewhat  arbitrarily  -  use  of  a  more  realistic  initial  covariance  matrix  as 
described  in  [3]  would  have  been  a  better  approach.  Additionally,  there  are  instances  where  the  filter 
loses  track  of  a  target,  for  example  for  Satellite  1  at  about  8  hours  into  the  run  there  are  4  time  steps 
for  which  the  solution  is  substantially  wrong  for  this  satellite.  Since  the  cardinality  stayed  at  5  targets, 
the  filter  essentially  was  tracking  a  clutter  point  during  this  time.  After  additional  measurements  were 
processed  the  filter  regained  the  correct  target. 
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Figure  12.  GEO  Element  Estimation  Errors  (Blue)  for  the  ASTRA  Cluster  with  3  Sigma 

Bounds  (Red) 
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5.0 


CONCLUSIONS 


Our  research  demonstrated  the  benefit  of  using  a  GEO  state  representation  in  support  of  SSA  for 
geostationary  space  objects.  We  found  that  linearized  EKF  implementations  of  both  the  ECI  and 
GEO  element  representations  were  suitable  for  GEO  angle  measurement  updates,  with  some  benefit 
to  the  GEO  elements  when  long  propagations  are  required.  We  used  both  Differential  Entropy  and 
the  Kullback-Liebler  measure  to  compare  the  four  filters. 

We  also  provided  an  initial  step  towards  the  application  of  a  GM-CPHD  filter  for  SSA  in  the  GEO 
regime.  The  approach  of  using  GEO  elements  was  shown  to  be  feasible  in  a  very  simple  simulation 
of  a  five  satellite  cluster,  observed  by  a  single  tracking  station,  with  significant  clutter  measurements. 
This  simulated  multi-target  scenario  is  a  good  starting  point  for  the  methods  developed  here. 
However,  there  is  substantial  room  for  improvement.  Future  efforts  must  consider  more  realistic 
initial  uncertainties,  measurement  intervals  and  duration,  better  models  of  the  measurement  errors, 
and  mis-modeling  of  the  spacecraft  motion.  After  the  end  of  this  project  period  we  were  able  to 
obtain  a  more  realistic  simulation  set  that  is  currently  being  processed  and  analyzed.  We  anticipate 
that  new  results  from  that  work  will  be  forthcoming. 

Given  the  promising  results  from  this  study,  we  are  eager  to  apply  the  methods  to  experimental 
observations  and  to  determine  if  the  GEO  approach  would  be  beneficial  in  reducing  the  computational 
requirements  for  large-scale  SSA  implementations. 
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LIST  OF  SYMBOLS,  ABBREVIATIONS,  AND  ACRONYMS 


CCAR 

Colorado  Center  for  Astrodynamics  Research 

CPHD 

Cardinalized  Probability  Hypothesis  Density 

DE 

Differential  Entropy 

ECEF 

Earth  Centered,  Earth  Fixed 

ECI 

Earth  Centered  Inertial 

EKF 

Extended  Kalman  Filter 

FISST 

Finite  Set  Statistics 

GEO 

Geostationary  Earth  Orbit 

GM 

Gaussian  Mixture 

GMM 

Gaussian  Mixture  Model 

GNSS 

Global  Navigation  Satellite  System 

JAT 

Java  Astrodynamics  Toolkit 

KL 

Kullback-Leibler 

MEO 

Medium  Earth  Orbit 

ODTBX 

Orbit  Determination  Toolbox 

OSPA 

Optimal  Sub-Pattern  Assignment 

PDF 

Probability  Density  Functions 

RIG 

Radial  In-track  Cross-track 

RK4 

Runge-Kutta  4th  Order  Integration  Scheme 

SRP 

Solar  Radiation  Pressure 

SSA 

Space  Situational  Awareness 

TLE 

Two-Line  Element 

UKF 

Unscented  Kalman  Filter 
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